Excitable Greenberg-Hastings cellular automaton model on scale-free networks 
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We study the excitable Greenberg-Hastings cellular automaton model on scale-free networks. We obtained 
analytical expressions for no external stimulus and the uncoupled case. It is found that the curves, the average 
activity F versus the external stimulus rate r, can be fitted by a Hill function, but not exactly, and there exists 
a relation F oc r a for the low-stimulus response, where Stevens-Hill exponent a ranges from a = 1 in the 
subcritical regime to a = 0.5 at criticality. At the critical point, the range reaches the maximal. We also 
calculate the average activity F k (r) and the dynamic range A fe (p) for nodes with given connectivity k. It is 
interesting that nodes with larger connectivity have larger optimal range, which could be applied in biological 
experiments to reveal the network topology. 
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Many systems in the real world, either naturally evolved 
or artificially designed, are organized in a network fashion 
[Q]] . The most studied networks are the exponential networks 
in which the nodes' connectivity distribution (the probabil- 
ity P(k) that a node if connected to other k nodes) is expo- 
nentially bounded. A typical model of this network is the 
Watts -Strogatz (WS) graph ||2|] which exhibits the "small- 
world"phenomenon that was observed in realistic networks 
On the contrary, it has been observed recently that the 
nodes' degree distribution of many networks have power-law 
tails (scale-free property, due to the absence of a characteristic 
value for the degrees) P(k) ~ k~~< with 2 < 7 < 3 fl. 

Scale-free (SF) networks have been widely studied in the 
past few years, mainly because of their connection to a lot of 
real-world structures, including networks in nature, such as 
the cell, metabolic networks and the food web, artificial net- 
works such as the Internet and the WWW, or even social net- 
works, such as sexual partnership networks [1]. The SF net- 
work is characterized with the power-law degree distribution. 
Thus there always exists a small number of nodes which are 
connected to a large number of other nodes in SF networks. 
This heterogeneity leads to intriguing properties of SF net- 
works. A lot of work has been devoted in the literature to 
the study of static properties of the networks, while many in- 
terests are growing for dynamical properties on these kind of 
networks. In the context of percolation, SF networks are sta- 
ble against random removal of nodes while they are fragile 
under intentional attacks targeting on nodes with high degree 
[5, 6]. Also, there are absences of epidemic thresholds in SF 
networks 10] and of the kinetic effects in reaction-diffusion 
processes taking place on SF networks JH]. 

In this paper, we will study the excitable Greenberg- 
Hastings cellular automaton (GHCA) model [9] on the SF 
network, especially we focus here on the Barabasi and Al- 
bert (BA) graph |4|. Due to experimental data which sug- 
gest that some classes of spiking neurons in the first layers 
of sensory systems are electrically coupled via gap junctions 
or ephaptic interactions, the GHCA model is employed to 
model the response of the sensory network to external stim- 
uli in some recent works. Two-dimensional deterministic cel- 
lular automaton model was studied by computer simulations 



ifToll . Analytical results have recently been obtained for the 
one-dimensional cellular automaton model under the two-site 
mean-field approximation [11]. In ref. fl2il . Kinouchi and 
Copelli studied the GHCA model on Erdos-Renyi random 
graph with stochastic activity propagation, and they found a 
new and important result that dynamic range maximized at 
the critical point, which is perhaps the first clear example of 
signal processing optimization at a phase transition. 

In the n-state GHCA model |9| for excitable systems, 
the instantaneous membrane potential of the i-th cell (i — 
1,...,N) at discrete time t is represented by Xi(t) £ 
{0, 1, . . . , n — 1}, n > 3. The state Xi(t) = denotes a neu- 
ron at its resting (polarized) potential, Xi(t) = 1 represents 
a spiking (depolarizing) neuron and Xi(t) — 2, . . . , n — 1 
account for the afterspike refractory period (hyperpolariza- 
tion). There are two ways for the i-th element to go from 
state Xi(t) = to Xi(t + 1) = 1: a) due to an external sig- 
nal, modelled here by a Poisson process with rate r (which 
implies a transition with probability A = 1 — exp(— r At) per 
time step); b) with probability p, due to a neighbor j being 
in the excited state in the previous time step. If Xi(t) > 1, 
then Xi(t + 1) = (xi(t) + l)mod n, regardless of the stimu- 
lus. In other words, the rules state that a neuron only spikes 
if stimulated, after which it undergoes an absolute refractory 
period before returning to rest. Time is discrete. We assume 
At = 1 ms which corresponds to the approximate duration of 
a spike and is the time scale adopted for the time step of the 
model. The number of states n therefore controls the duration 
of the refractory period (which corresponds to n — 2, in ms). 
In the biological context, r could be related for example with 
the concentration of a given odorant presented to an olfactory 
epithelium lfl3ll . or the light intensity stimulating a retina (3. 
We shall refer to r as the stimulus rate or intensity. 

BA graph is a kind of SF networks and can be constructed 
according to ref. 0]. Starting from a small number too 
of nodes, every time step a new vertex is added, with m 
links that are connected to an old node i with probability 
= ki/^2j kj, where ki is the connectivity of the ith node. 
After iterating this scheme a sufficient number of times, we 
obtain a network composed by N nodes with connectivity dis- 
tribution P(k) ~ fc -3 and average connectivity (k) = 2m. 
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For this kind network, the absence of a characteristic scale for 
the connectivity makes highly connected nodes statistically 
significant, and induces strong fluctuations in the connectiv- 
ity distribution which cannot be neglected. 

Let pt(s) be the densities of neurons which are in state s at 
time t. We have the normalization condition, Y^"=o P*( s ) ~ 
1. Since the dynamics of the refractory state is deterministic, 
the equations for s > 2 are simply 
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Pt+i(n - 1) = p t (n - 2) . 
By imposing the stationarity condition, we have 

p t (0) = l-(n-l)p 4 (l), 



(1) 



(2) 



Following ideas developed by Pastor-Satorras et al. |7|], to 
take into account strong fluctuations in the connectivity distri- 
bution, we consider the relative density ( o f fe (l) of active nodes 
with given connectivity k by the following equation 

d tPt k {\) = -p t k (l) + X Pt k (0) + kp(l - X) Pt k (0)Q( Pt (1)), 

(3) 

where Q(p t (l)) is the probability that any given link point 
to an active node and is assumed as a function of the total 
density of exciting nodes pt(l). In the steady state, pt(l) is 
just a function of A and p. Thus, the probability 6 becomes 
an implicit function of A and p. By imposing the stationarity 
condition, 9 t p t fe (l) = 0, we obtain 

Jb m= A+(l-A)pfc9(A,p) 

9 1 ' 1 + (n - 1)[A + (1 - X) P kQ(X,p)} ' k ' 

This set of equations show that the higher the node connec- 
tivity, the higher the probability to be in a spiking state. This 
inhomogeneity must be taken into account in the computa- 
tion of Q(X,p). Indeed, the probability that a link points to 
a node with q links is proportional to qP(q). In other words, 
a randomly chosen link is more likely to be connected to an 
exciting node with high connectivity, yielding the relation 



e(A, P ) = ]T 



kP(k) P k (i) 



(5) 



Since p k (l) is on its turn a function of 0(A,p), we obtain a 
self-consistency equation that allows to find Q(X,p) and an 
explicit form for Eq. (@}. Finally, we can evaluate the order 
parameter (persistence) p(l) using the relation 



p(l) = ^P(fc)p fc (l). 

k 



(6) 



For the BA graph, the full connectivity distribution is given 
by P(k) — 2m 2 k~ 3 , where m is the minimum number of 
connection at each node. By noticing that the average con- 
nectivity is (k) = kP(k)dk = 2m, Eq. © gives 



&(X,p) = m 



X+{l-X)pkQ(X,p) 



k 2 1 + (n - 1)[A + (1 - X)pke{X,p)} ' 

(7) 
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FIG. 1 : In the absence of stimulus, the average activity Fasa func- 
tion of 1/p for BA networks of size TV = 10 4 and mo = m = 4. The 
linear behavior on the semi-logarithmic scale proves an exponential 
behavior predicted by Eq. ( |10| l. 



which yields the solution 

e(A lP ) = £ + d - - 1 
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maO(X,p) + b' 
(8) 

where a — (n — 1)(1 — X)p and b = 1 + (n — 1)A. We 
can solve the above equation to obtain the solution Q(X,p). 
Theoretically, combing Eqs. (0]l , © and ([§), we can obtain 
the final result of p(l). In some special cases the analytical 
expressions can be obtained. 

i) no external stimulus, i.e., A = or r = 0. We can obtain 
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(9) 



Combing Eqs. ©, ©, and (O, we find the solution of the 
density of active nodes when there is no external stimulus, 



P(l) 



n-l 



- 1/ rap 



(10) 



ii) the uncoupled case, i.e., p — 0. Combing Eqs. and 
©, we have 



l + ( n -l)A' 



(11) 



which is independent of the network's topology, so it is also 
have been obtained in other networks, such as one dimen- 
sional case [11]. When the external stimulus intensity r is 
very small, p(l) ~ r _1 . 

We define the average activity 
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FIG. 2: Response curves (mean firing rate vs. stimulus rate). Points 
represent simulation results with TV = 10 4 , mo = m — 4, n — 5 
states and T = 10 3 ms, from p — to p — 0.12 (in intervals of 
0.02). These curves are power laws F oc r a with a = 1 (subcritical) 
and a = 1/2 (critical). Inset: spontaneous activity Fo vs. branching 
probability p, the critical point is p c = 0.06. 



where T is a large time window (of the order of 10 4 time 
steps). In the stationary state, it is obviously that F = p(l). 
To confirm the picture extracted from the above analytic treat- 
ment, we perform numerical simulations on the BA network. 
Figure Q] shows the behavior of the average activity F in the 
absence of stimulus. All the plots decays with an exponent 
form F ~ exp(— c/mp), where c is a constant. The numer- 
ical value obtained c = 0.226 is in good agreement with the 
theoretical prediction c = 1/m = 0.25. 

In Fig. |2] we show the average activity F versus the exter- 
nal stimulus rate r for different branching probability p. The 
inset shows the spontaneous activity Fo versus the branching 
probability p in the absence of stimulus (r — 0), there is a crit- 
ical point p c = 0.06, only at p > p c there is a self-sustained 
activity, F > 0. The curves F(r) could be fitted by a Hill 
function, but are not exactly Hill, and there exists the rela- 
tion F oc r a for the low-stimulus response. It is found that 
the Stevens-Hill exponent a changes from a = 1 in the sub- 
critical regime to a — 0.5 at criticality. This important point 
was also reported by ref. [ 12] where the network is an Erdos- 
Renyi random graph. We also notice that apparent exponents 
between 0.5 and 1.0 are observed ifTTll if finite size effects are 
present, that is, if N is small. 

As a function of the stimulus intensity r, networks have 
a minimum response Fo (= for the subcritical and critical 
cases) and a maximum response F max (due to the absolute 
nature of the refractory period, F max = 1/n, which can be 
obtained by setting A = 1 in Eqs. (fTll). Ref. [ 12] defines the 
dynamic range A = 10 log(ro.g/ro.i) as the stimulus interval 
(measured in dB) where variations in r can be robustly coded 
by variations in F, discarding stimuli which are too weak to 
be distinguished from Fo or too close to saturation. The range 
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FIG. 3: Dynamic range vs. branching probability p. The curve is ob- 
tained by calculating the data from Fig(2] In the subcritical regime, 
the dynamic range A(p) increases monotonically with p, In the su- 
percritical regime, A(p) decreases when p increases. There is a max- 
imal range precisely at the critical point p c = 0.06. 



[ r o. i, r o.9] is found from its corresponding response interval 
[Fo.i,F . 9 ], where F x = F + x(F max - F ). 

Figure [3] depicts the dynamic range A versus the branch- 
ing probability p. In the subcritical regime, the dynamic 
range A(p) increases monotonically with p, In the supercrit- 
ical regime, A(p) decreases when p increases. There is a 
maximal range precisely at the critical point. This result was 
also found in Erdos-Renyi random graph 11211 . Kinouchi and 
Copelli explained the result as "In the subcritical regime, sen- 
sitivity is enlarged because weak stimuli are amplified due to 
activity propagation among neighbours. As a result, the dy- 
namic range A(p) increases monotonically with p. In the 
supercritical regime, the spontaneous activity Fo masks the 
presence of weak stimuli, therefore A(p) decreases. The opti- 
mal regime occurs precisely at the critical point" lfl2ll . We also 
perform simulations on other complex networks, such as WS 
networks and random SF networks with different 7, and ob- 
tain the same behavior of the average activity F as it is shown 
in Fig. [2] the optimal regime occuring precisely at the critical 
point. We state that this phenomenon is the universal behav- 
ior of excitable GHCA model on complex networks and the 
explanation provided by Kinouchi and Copelli is also suitable 
for other types of complex networks. 

In SF networks, there exist strong fluctuations in the con- 
nectivity distribution. It is worthy to investigate the behavior 
of the average activity F k for nodes with given connectivity 
k. In Fig. [4]we plot the quantity F k versus the external stimu- 
lus intensity r for different branching probability p. Figs. Ufa) 
and (b) correspond to the node's connectivity k = 4 and 36, 
respectively. It is found that the curves F k (r) have the similar 
property of F(r) which was shown in Fig. |2l i.e., the Stevens- 
Hill exponent a changes from a = 1 in the subcritical regime 
to a = 0.5 at criticality. 

Finally, we calculate the dynamic range A k (p) of nodes 
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FIG. 4: Average activity F k for nodes with given connectivity k vs. 
external stimulus rate r. (a) for connectivity k — 4, and (b) for 
connectivity k — 36. Simulations are performed with N = 10 4 , 
mo = m = 4, n = 5 states and T = 10 4 ms, from p = to p = 
0.12 (in intervals of 0.02). These curves are power laws F k oc r a 
with a = 1 (subcritical) and a = 1/2 (critical). 
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with given connectivity k, and the result is shown in Fig. [5] 
The phenomenon that the optimal regime occurs precisely at 
the critical point recurs. It is notable that the optimal dynamic 
range for nodes with given connectivity k increases with k, 
i.e., for the node with larger connectivity, its corresponding 
larger optimal range. One can investigate every node's opti- 
mal dynamic range and calculate their fluctuations, since the 
fluctuations of node's optimal dynamic range reflect the fluc- 
tuations of the node's connectivity in the network. To some 
extent, we can discovery the topology of the network by in- 
vestigating the dynamics on it. 

In summary, we have investigated the behavior of the ex- 
citable GHCA model |9] on BA networks. We found out that 
the curves of the average activity F as a function of the ex- 
ternal stimulus rate r can be fitted by a Hill function, but are 
not exactly Hill, and there exists a relation F oc r a for the 
low-stimulus response. The Stevens-Hill exponent a changes 
from a — 1 in the subcritical regime to a — 0.5 at criticality. 
There is a maximal range precisely at the critical point. We 
also observed these results numerically in other kind of com- 
plex networks. We conclude that these phenomena are the 
universal behaviors of the excitable GHCA model on com- 
plex networks. Due to strong fluctuations in the connectivity 
distribution on the BA graph, we calculated the average activ- 
ity F k (r) and the dynamic range A k (p) for nodes with given 
connectivity k. The two quantities F k (r) and A k (p) have the 
similar behavior as that of F(r) and A(p), respectively. It 
is interesting that nodes with larger connectivity have larger 
optimal range. This property could be applied in biological 
experiment revealing the network topology. 
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FIG. 5: Dynamic range A fc (p) of nodes with given connectivity k 
vs. branching probability p. Points represent the results calculated 
from simulation with N = 10 4 sites, mo = m = 4, n = 5 states 
and T = 10 4 ms. Dynamic range A k (p) is optimized at the critical 
point p c = 0.06, although the connectivity k is different . 
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